!
!     Kombine     Combine the set of frames.
!     Copyright (C) 2000  Filip Hroch, Masaryk University, Brno, CZ
!
! $Id$
!

       Program ProgKombine
!
! Experimental version 
!
! Sum (combine) the frames.
! Using the output of the munimatch (with -t switch), eg. 'TRA'
! files with format:
! image1 image2 scale dx1, dx2 cos(fi) sin(fi) sin(fi) cos(fi)

     implicit none

! Parameters:

!    Version:
         Character(len=80),parameter :: version = 'Kombine v. 0.0.3'

! Initial parameters.
!    Output image name:
       Character(len=256) :: namekombi='kombine.fits'

!    Minimum of values
       Real :: minvalue = 0.0

!    Maximum of values
!       Real :: maxvalue = 65535.0

!-------

       Character(len=666) :: sb, sbtra
       Character(len=80) :: buf

       Logical :: simple,extend,anyf,exist
       Integer :: n, nalldate, maxalldate
       Integer :: i, j, istat, ostat ,bitpix, naxis, naxes(2) 
       Integer :: BlockSize,pcount,gcount

! image dimensions
       Integer ::  nx = -1, ny = -1 
! image
       Real, allocatable :: KOMBINE(:,:), IMAGE(:,:)
       Real :: bscale, bzero
       Real :: oneexp, exptime
       character(len=80) :: onedate,object
       character(len=80),allocatable :: alldate(:), alld(:)
       Real :: xx,yy,c,r(2),m(2,2),t1(2),t2(2)
       integer :: ix,jy

!------

!       write(*,*) trim(version)

! init sum variables
       object = ''
       exptime = 0.0
       nalldate = 0
       maxalldate = 2
       allocate(alldate(maxalldate))

       ! input of top line with init values
       read(*,*,end=666,iostat=j) sb
       if( j /= 0 ) stop 'Bad format of input data.'
       if( sb /= '' ) namekombi = sb

       istat = 0

       n = 0
! nacteni jmena souboru
       Do
         Do
           Read(*,*,END=20,ERR=666) sb
           if( sb /= ' ' ) exit              ! blank line
         enddo
         write(*,*) trim(sb),':'
         
         istat = 0
! input initialisation of the current file
         Call ftopen(25,sb,0,BlockSize,istat)
         Call ftghpr(25,2,simple,bitpix,naxis,naxes,pcount,gcount,extend,istat)
         if( istat /= 0 ) goto 666
         CALL ftgkye(25,'BSCALE',bscale,buf,istat)
         CALL ftgkye(25,'BZERO',bzero,buf,istat)
!        BSCALE, BZERO not presented
         if( istat /= 0 ) then
            istat = 0
            bscale = 1.0
            bzero = 0.0
         endif
         call ftgkye(25,'EXPTIME',oneexp, buf, istat)
         if( istat /= 0 ) then
            oneexp = -1.0
            istat = 0
         endif
         call ftgkys(25,'DATE-OBS',onedate, buf, istat)
         if( istat /= 0 ) then
            onedate = ''
            istat = 0
         endif
         call ftgkys(25,'OBJECT',object, buf, istat)
         if( istat /= 0 ) then
            object = ''
            istat = 0
         endif

! first image - internal buffers initialisation
         if( nx < 0 .and. ny < 0 )then
            nx = naxes(1)
            ny = naxes(2)
            Allocate( KOMBINE(nx,ny) )
            Allocate(   IMAGE(nx,ny) ) 
            IMAGE = 0.0
            KOMBINE = 0.0
         else       ! checking image dimensions
            if( nx /= naxes(1) .or. ny /= naxes(2) ) &
              stop 'Different dimensions of different images.'
         endif 

! image - data input
         n = n + 1
         call FTG2DE(25,1,minvalue,nx,nx,ny,IMAGE,anyf,istat)
         if( istat /= 0 ) goto 666
         call ftclos(25,istat)
         if( istat /= 0 ) goto 666

! image - transformation input
         sbtra = sb(:index(sb,'.f')+index(sb,'.F'))//'TRA'
         open(27,file=sbtra,status='old',iostat=i)
         if( i == 0 ) read(27,*,iostat=i) sb,sb,c,r,m,t1,t2
         if( i /= 0 ) then
            write(*,*) 'No such TRA file or invalid format. Skipped.'
         endif
         close(27)
         if( i == 0 ) then

            ! compute sumation
            do i = 1, nx
               do j = 1, ny
                  xx = i - t1(1)
                  yy = j - t1(2)
                  ix = nint(c*(r(1) + xx*m(1,1) + yy*m(1,2)) + t2(1))
                  jy = nint(c*(r(2) + xx*m(2,1) + yy*m(2,2)) + t2(2))
                  if( 1<=ix .and. ix<=nx .and. 1<=jy .and. jy<=ny ) then
                     KOMBINE(i,j) = KOMBINE(i,j) + IMAGE(ix,jy)
                  endif
               enddo
            enddo

            if( onedate /= '' .and. oneexp > 0.0 ) then
               write(buf,*) trim(onedate),oneexp
               exptime = exptime + oneexp

               nalldate = nalldate + 1
               if( nalldate > maxalldate ) then       ! need to write after tail
                  allocate(alld(maxalldate))
                  alld = alldate
                  deallocate(alldate)
                  maxalldate = 2*maxalldate
                  allocate(alldate(maxalldate))
                  alldate(1:maxalldate/2) = alld
                  deallocate(alld)
               endif
               alldate(nalldate) = buf 

            endif
         end if
      enddo
 20   continue

      if( n == 0 ) stop 'No input image.'
      
! form the output image
       ostat = 0
       Inquire(file=namekombi,exist=exist)
       if( exist ) then
         open(26,file=namekombi)
         close(26,status='DELETE')
       endif
       Call ftinit(26,namekombi,1,ostat)
       naxis = 2
       naxes = (/nx, ny/)
       bitpix = 32
       Call ftphps(26,bitpix,naxis,naxes,ostat)
       Call ftpkys(26,'OBJECT',object,'',ostat)
!       Call ftpkye(26,'BSCALE',bscale,10,'',ostat)
!       Call ftpkye(26,'BZERO',bzero,10,'',ostat)
       call ftpkye(26,'EXPTIME',exptime,10,'Sum of all times',ostat)    
       if( nalldate > 0 ) then
          write(buf,*) 'The image is a result of composition of ',nalldate,' exposure(s).'
          call ftpcom(26,buf,ostat) 
          call ftpcom(26,'The time of start and the exposure time for each image used:',ostat)
          do i = 1, nalldate
            call ftpcom(26,alldate(i),ostat)
          enddo
          call ftpcom(26,'All times are derived from DATE-OBS and EXPTIME keywords.',ostat)
       endif
       if( ostat /= 0 ) goto 666

       Call ftp2de(26,1,nx,nx,ny,KOMBINE,ostat)

       CALL ftclos(26,ostat)
       if( ostat /= 0 ) goto 666

! print output info
       write(*,*) 'Output image: ',trim(namekombi)
       write(*,*) 'Dimension: ',nx,'x',ny
!       write(*,*) 'Bye.'
       stop 'Bye.'

666    continue 

       if( istat +  ostat /= 0 )then
        write(*,*) 'The following error(s) has been occured:'
        do
          Call ftgmsg(buf)
          if( buf == ' ' ) exit
          write(*,*) trim(buf)
        enddo
      endif

      stop 'Crached...'

      End       

